{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge Sort\n",
    "\n",
    "Code examples from [Think Complexity, 2nd edition](http://greenteapress.com/wp/complexity2), Appendix A\n",
    "\n",
    "Copyright 2017 Allen Downey, [MIT License](http://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import string\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import thinkplot\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical order of growth\n",
    "\n",
    "Sometimes we can figure out what order of growth a function belongs to by running it with a range of problem sizes and measuring the run time.\n",
    "\n",
    "To measure runtimes, we'll use `etime`, which uses `os.times` to compute the total time used by a process, including \"user time\" and \"system time\".  User time is time spent running your code; system time is time spent running operating system code on your behalf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def etime():\n",
    "    \"\"\"Measures user and system time this process has used.\n",
    "\n",
    "    Returns the sum of user and system time.\"\"\"\n",
    "    user, sys, chuser, chsys, real = os.times()\n",
    "    return user+sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`time_func` takes a function object and a problem size, `n`, runs the function, and returns the elapsed time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def time_func(func, n):\n",
    "    \"\"\"Run a function and return the elapsed time.\n",
    "    \n",
    "    func: function\n",
    "    n: problem size\n",
    "    \n",
    "    returns: user+sys time in seconds\n",
    "    \"\"\"\n",
    "    start = etime()\n",
    "    func(n)\n",
    "    end = etime()\n",
    "    elapsed = end - start\n",
    "    return elapsed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`run_timing_test` takes a function, runs it with a range of problem sizes, and returns two lists: problem sizes and times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_timing_test(func, max_time=1):\n",
    "    \"\"\"Tests the given function with a range of values for n.\n",
    "    \n",
    "    func: function object\n",
    "\n",
    "    returns: list of ns and a list of run times.\n",
    "    \"\"\"\n",
    "    ns = []\n",
    "    ts = []\n",
    "    for i in range(6, 28):\n",
    "        n = 2**i\n",
    "        t = time_func(func, n)\n",
    "        print(n, t)\n",
    "        if t > 0:\n",
    "            ns.append(n)\n",
    "            ts.append(t)\n",
    "        if t > max_time:\n",
    "            break\n",
    "\n",
    "    return ns, ts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`fit` takes the lists of ns and ts and fits it with a curve of the form `a * n**exp`, where `exp` is a given exponent and `a` is chosen so that the line goes through a particular point in the sequence, usually the last. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fit(ns, ts, exp=1.0, index=-1):\n",
    "    \"\"\"Fits a curve with the given exponent.\n",
    "    \n",
    "    ns: sequence of problem sizes\n",
    "    ts: sequence of times\n",
    "    exp: exponent of the fitted curve\n",
    "    index: index of the element the fitted line should go through\n",
    "    \n",
    "    returns: sequence of fitted times\n",
    "\n",
    "    \n",
    "    \"\"\"\n",
    "    # Use the element with the given index as a reference point, \n",
    "    # and scale all other points accordingly.\n",
    "    nref = ns[index]\n",
    "    tref = ts[index]\n",
    "\n",
    "    tfit = []\n",
    "    for n in ns:\n",
    "        ratio = n / nref\n",
    "        t = ratio**exp * tref\n",
    "        tfit.append(t)\n",
    "\n",
    "    return tfit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`plot_timing_test` plots the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_timing_test(ns, ts, label='', color='blue', exp=1.0, scale='log'):\n",
    "    \"\"\"Plots data and a fitted curve.\n",
    "\n",
    "    ns: sequence of n (problem size)\n",
    "    ts: sequence of t (run time)\n",
    "    label: string label for the data curve\n",
    "    color: string color for the data curve\n",
    "    exp: exponent (slope) for the fitted curve\n",
    "    \"\"\"\n",
    "    tfit = fit(ns, ts, exp)\n",
    "    plt.plot(ns, tfit, color='0.7', linewidth=2, linestyle='dashed')\n",
    "    plt.plot(ns, ts, 's-', label=label, color=color, alpha=0.5, linewidth=3)\n",
    "    plt.xlabel('Problem size (n)')\n",
    "    plt.ylabel('Runtime (seconds)')\n",
    "    plt.xscale(scale)\n",
    "    plt.yscale(scale)\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For small values of `n`, the runtime is so short that we're probably not getting an accurate measurement of just the operation we're interested in.  But as `n` increases, runtime seems to converge to a line with slope 1.  \n",
    "\n",
    "That suggests that performing append `n` times is linear, which suggests that a single append is constant time.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing sort algorithms\n",
    "\n",
    "NumPy provides implementations of three sorting algorithms, quicksort, mergesort, and heapsort.\n",
    "\n",
    "Read about each of these algorithms to see what order of growth they belong to.\n",
    "\n",
    "Now let's see if we can characterize their asymptotic behavior.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def test_quicksort(n):\n",
    "    xs = np.random.normal(size=n)\n",
    "    xs.sort(kind='quicksort')\n",
    "\n",
    "ns, ts = run_timing_test(test_quicksort)\n",
    "plot_timing_test(ns, ts, 'test_quicksort', exp=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quicksort is hard to distinguish from linear, up to about 10 million elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def test_mergesort(n):\n",
    "    xs = np.random.normal(size=n)\n",
    "    xs.sort(kind='mergesort')\n",
    "\n",
    "ns, ts = run_timing_test(test_mergesort)\n",
    "plot_timing_test(ns, ts, 'test_mergesort', exp=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merge sort is similar, maybe with some upward curvature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def test_heapsort(n):\n",
    "    xs = np.random.normal(size=n)\n",
    "    xs.sort(kind='heapsort')\n",
    "\n",
    "ns, ts = run_timing_test(test_quicksort)\n",
    "plot_timing_test(ns, ts, 'test_heapsort', exp=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three methods are effectively linear over this range of problem sizes.\n",
    "\n",
    "And their run times are about the same, with quicksort being the fastest, despite being the one with the worst asympotic performance in the worst case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing Merge Sort\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def merge_sort_norec(xs):\n",
    "    N = len(xs)\n",
    "    left = xs[:N//2]\n",
    "    right = xs[N//2:]\n",
    "    \n",
    "    left.sort()\n",
    "    right.sort()\n",
    "    \n",
    "    return merge(left, right)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This version breaks the array in half, uses `np.sort` to sort the two halves, then uses merge to put the halves together.\n",
    "\n",
    "**Exercise:** Write a function called `merge` that takes two sorted NumPy arrays, `left` and `right`, and returns a new array that contains all elements from `left` and `right`, sorted.  (where \"sorted\" means in ascending order, or non-decreasing, to be more precise).\n",
    "\n",
    "Note: this function is not hard to write, but it is notoriously difficult to get all of the edge cases right without making the function unreadable.  Take it as a challenge to write a version that is correct, concise, and readable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xs = np.random.random(10)\n",
    "ys = np.random.random(10)\n",
    "xs.sort()\n",
    "ys.sort()\n",
    "res = merge(xs, ys)\n",
    "all(sorted(res) == res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Starting with `merge_sort_norec`, write a function called `merge_sort_rec` that's fully recursive; that is, instead of using `numpy.sort` to compute the DFTs of the halves, it should use `merge_sort_rec`.  Of course, you will need a base case to avoid an infinite recursion.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your method by running the code in the next cell, then use `test_merge_sort_rec`, below, to check the performance of your function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xs = np.random.random(10)\n",
    "\n",
    "res = merge_sort_rec(xs)\n",
    "all(sorted(res) == res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def test_merge_sort_rec(n):\n",
    "    xs = np.random.normal(size=n)\n",
    "    spectrum = merge_sort_rec(xs)\n",
    "\n",
    "ns, ts = run_timing_test(test_merge_sort_rec)\n",
    "plot_timing_test(ns, ts, 'test_merge_sort_rec', exp=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If things go according to plan, your implementation of merge sort should be close to linear, or a little steeper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
